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Abstract 

A simplified numerical model of the plasma discharge within the SPT- 100 stationary plasma thruster 
was developed to aid in understanding thruster operation. A one dimensional description was used. 
Non-axial velocities were neglected except for the azimuthal electron velocity. A nominal operating 
condition of 4.5 mg/s of xenon anode flow was considered with 4.5 Amperes of discharge current, and a 
peak radial magnetic field strength of 130 Gauss. For these conditions, the calculated results indicated 
ionization fractions of 0.99 near the thruster exit with a potential drop across the discharge of 
approximately 250 Volts. Peak calculated electron temperatures were found to be sensitive to the choice 
of total ionization cross section for ionization of atomic xenon by electron bombardment and ranged from 
51 eV to 60 eV. The calculated ionization fraction, potential drop, and electron number density agree 
favorably with previous experiments. Calculated electron temperatures are higher than previously 
measured. 


Nomenclature 


A Cross sectional area of discharge chamber 
A Coefficient matrix 

ao Bohr radius 

— £ 

b Column vector of inhomogeneous terms 

B r Radial magnetic field strength 

D e Electron diffusion coefficient 
£f Ionization energy 

Ionization energy of hydrogen 

— > 

E Electric field 

e Elementary charge 

f(E) Electron energy distribution function 
J Discharge current 

Je Electron current density 

k Boltzmann's constant 

fcj Volumetric rate of ionization 
L Discharge length (anode to cathode) 

m Mass flow rate 

m a Atom mass 

m e Electron mass 

mi Ion mass 

ria Atom number density 

tie Electron number density 

«i Ion number density 


P e Electron pressure 
Q(E) Collision cross section 
q e Conductive heat flux due to electrons 
T a Atom temperature 
T e Electron temperature 

T,- Ion temperature 
— £ 

U e Electron drift velocity 
U e Axial electron drift velocity 
lI e Q Circumferential electron drift velocity 
U a Electron drift velocity 
Ui Ion drift velocity 

w Net rate of electron/ ion production 

z Spatial coordinate 

x Column vector of unknowns 

£i Ionization energy 

He Electron mobility 

< o Average collision cross section 
v«j Atom ion momentum transfer collision 
frequency 

v ea Electron atom momentum transfer collision 
frequency 

v e i Electron ion momentum transfer collision 
frequency 

<p Electric potential 

iction 


During stationary plasma thruster operation a plasma is sustained within an annular discharge chamber by 
an axial electric field established between an external hollow cathode and an anode located at the rear. 
Near the anode, electron-ion pairs are formed through collisional ionization via electron bombardment. 
The interaction between the electrons and the crossed axial electric field and radial magnetic field 
established by electromagnets redirects the electrons in a azimuthal direction greatly reducing the electron 
conductivity in the axial direction. The reduced axial electron conductivity permits large electric fields 
which efficiently accelerate thrust producing ions. A cross sectional schematic of the discharge chamber 
is shown in Figure 1. The objective of this investigation was to develop a tool capable of examining the 
details of these processes using the minimum subset of the governing equations required to retain all 
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essential physical phenomenon. As a result of this 
approach, non-axial velocities were neglected except for 
azimuthal electron velocity. These limitations prevent 
direct insight into important issues such as thruster lifetime 
and spacecraft integration associated with non-axial ion 
velocities but allow numerical solutions of the discharge 
physics to be obtained using a personal computer in a 
timely fashion. 

Specifically, the model consisted of seven coupled 
algebraic and differential equations. These included ion 
continuity, global continuity, current continuity, ion 
momentum, electron momentum and electron energy. 
Assumptions included one dimensional flow, quasi- 
neutrality, no doubly charged ions, constant and equal 
atom and ion temperatures, and constant atom velocity. 
Seven unknowns were calculated: electron/ion density, 
atom density, ion velocity, electron velocity, electron 
temperature, and potential. Initial conditions included 
discharge current, mass flow rate, initial ionization fraction 
and initial electron temperature. The axial distribution of 
the radial magnetic field was also assumed. This approach 
differs from that of past investigators. Melikov'sfl] one- 
dimensional model relied on experimentally determined 
values for the potential and electron temperature 
distribution. Lentz [2] relied on Bohm diffusion to 
describe the axial motion of the electrons across the 
magnetic field lines. More recently, Fife[3] has 
successfully extended this work to include both axial and 
radial dimensions. 



Figure 2; Cross Sectional Schematic of the Annular 
Section of the SPT Discharge Chamber and 
Assumed Radial Magnetic Field Distribution. 


In this paper, model formulation is presented including additional assumptions and simplifications not 
previously mentioned. The first-order explicit solution technique is also presented. Calculated results are 
presented as is the sensitivity of these results to the initial conditions and inputted parameters. Finally, 
where appropriate, the results are compared to experimental data. 


Analysis 

General forms of the energy, momentum, and continuity equations served as the starting point for this 
model. These general equations are presented. Assumptions and simplifications required to obtain the 
one dimensional scalar equations in the axial direction which were solved are stated. The one dimensional 
scalar equations are also presented. In all cases a singly ionized, steady-state, quasi-neutral plasma was 
considered (i.e. n e =ni). Neglecting radiation, viscous disipation, and flow work the general form of the 
electron energy equation considered was[4]: 


”“■> riigUg — & 5 — > — > 2 itu o 

V • ( U e n e 5 + Ue n e T k T e ) = - V.q e + J e • E - £ — - v e ; ii*. ~ k(T e -Tj) - w la) 

j 

The convective flux of kinetic energy on the right hand side of the previous equation included the flux of 

azimuthal kinetic electron energy (i.e. U^= v^ z +u^). Therefore, neglecting heat conduction, the electron 
energy equation was written as: 
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I ne me <l£ + i +U ^ + f He Hek-gf - He e U e ^ = 

^ea^ekCTg-Ta ) - 3 v ei n ek(T e -Tj)-(£j+2kT e ) w - ngmgUgUey 


(lb) 


This was the form of the electron energy equation which was solved. Similarly, in general form, 
neglecting the ion pressure gradient, electron-ion and atom-ion collisional momentum exchange] 
gravitational forces, shear stresses, and magnetization of ions, the ion momentum equation was given by:* 


nemiUi*VUi = eniE - irtiUiW 


or in one dimension: 




d Ui 


ne e + nemjUi-^- = - m;Uiw 


(2a) 


(2b) 


The previous assumptions were made on the basis that the dominate force on the ions is the electric field 
component of the Lorentz force. Identically, without simplification, the steady state ion continuity was: 


or in one dimension: 


v. (nj Uj )= w 

„ iUi tt dn e 
n “dz +u >^ - w 


( 3 a) 


(3b) 


The electron momentum equation was considered neglecting shear stresses and gravitational forces. In 
this case the general form of the electron momentum vector equation was[4]: 


— ) — ) 


rig mg ( U e • V U e ) + VP e = -eneE + J e x B -XneitieVejdJe-Uj)- mgUeW 

or in one dimension: 

. dT e dne 4Ug dip 

^ it + kTe dT + ^Ue-sr ' ^ e dl = 

nedJee B r - nemeVea(U e -Ua) - nemeVei(U e -Ui) - meU e w 


(4a) 


(4b) 


Again, although this equation is the axial component of the electron momentum equation, it includes the 

azimuthal electron velocity. This term, resulting from the J e x B Lorentz term is thought to be critical 
for accurately modeling the axial motion of electrons. However, the radial magnetic field and azimuthal 
electron velocities are treated as input parameters to the model. 

This differs from an alternate approach [2,3] utilized by previous investigators who neglect the Lorentz 
force due to the magnetic field as well as the convective term and the electron production term. For this 
case, the electron momentum equation became: 


VP e = -eneE + 1^111^ (U e -Ui) 


(5) 


Assuming the ion drift velocity is negligible with respect to the electron drift velocity, the electron 
mobility and a electron diffusion coefficient were defiened such that: 


VP P 


where, 


He — ■ Me E - D e p 


Me = 


and D e = 


e 

kT P 


(6) 


(7,8) 


m^Vgi c xTieV e j 

For cases where there are strong magnetic fields, the classical diffusion coefficient has been shown to be 
insufficient to correctly describe electron motion. In this case, the empirical Bohm diffusion coefficient 
was used. 
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D, 


Bohm ! 


kT e 
16e B 


(9) 


This approach relaxes the necessity of knowing the azimuthal electron velocity, decoupling the axial 
component of the electron momentum equation from the azimuthal component. However, the solution 
relies on an empirical diffusion coefficient that has not been shown to be applicable in this case. 

The last differential equation used in this model to form a solution set was the atom momentum equation. 
The generalized form used neglected gravitational forces, shear stresses, atom pressure, electron-atom 
collisional momentum exchange, and the net rate of production of atom momentum due to ionization and 
recombination : 


or in one dimension: 


n a m a (U a • V U a ) = - n a m a Vai (U a -Ui) 


d U a 

n a m a U a -^~ = - n a m a Vai(U a -Ui) 


(10a) 


(10b) 


The term accounting for atom-ion mometum exchange, which was retained in equation 10, was neglected 
in the ion momentum equation because the amount of momentum transfered to the ions by collisions with 
atoms is negligible with respect to the acceleration of the ions by the electric field. Two remaining 
algebraic equations were used to form a closed set of equations: current continuity and global continuity. 


% eUiA - rig eU e A = J 


( 11 ) 


m 


iWiTaUa +nemjUi = j 


( 12 ) 


For this investigation A was considered constant throughout the computational domain, which does not 
accurately descide the expansion takes place downstream of the discharge chamber exit. Equations lb, 
2b, 3b, 4b, 10b, 11, and 12 form a set of seven coupled algebraic and differential equations with seven 

unknowns: U e , Ui, U a , T e , n e , n a , <t>. Magnetic field strength and azimuthal electron velocity appear as a 
product in the relevant equations, and as a result were treated as a single lumped parameter. The net rate 
of production of electrons and ions and the collisional momentum exchange frequencies were calculated 
in the following fashion, w, the net rate of electron/ion production is calculated based on the rate of 
ionization from electron impact and the rate of recombination due to collisions with the discharge chamber 
walls (based on a Bohm velocity). 


w = 


ll/2 


(13) 


Because the computational domain extends beyond the exit of the discharge chamber, the loss term is set 
to zero outside the discharge chamber. The rate of ionization, is calculated based on the following 
integral: 


ki = J VE f(E) Q(E)dE 


(14) 


where f(E) was assumed to be a Maxwellian electron energy distribution function, i.e.: 

((E) = q= (j^) 3/2 E 1/2 exp( -E/kT e ) (15) 

and Q(E), the collision cross section for ionization, was initially based on the empirical formulation 
presented by Drawin[5]. 


4 



Q(E) = 2.6671 ag 


.'jz ln(1.25E/Ej) dE 


(16) 


Experimentally measured cross section were alternatively used[6,7]. The numerically integrated rate 
constants obtained using equation (14) are shown in Figure 2 based on the empirical collision cross 
sections and the experimental values. While the experimental cross sections seem significantly different, 
the calculated rates show a much smaller difference. The effect of using these different rates will be 
discussed along with the results. 
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Figure 2: Total ionization cross section for ionization of atomic xenon by electron bombardment and 
calculated rates of ionization assuming an Maxwellian electron energy distribution function. 


The electron-ion momentum transfer collision frequency was taken from Mitchner and Kruger[8]. This 
expression was: 

For electron-atom and atom-ion momentum transfer collisions the geometric cross section was used. 


( 8kT e \l/2 , 

) and v a i = nj <o> Uj 

70IV, / 


(18,19) 


In all cases charge exchange collisions within the discharge were neglected. The previously mentioned 
closed set of equations was solved in the following fashion. was eliminated from eqns. 1-4 using 

eqn. 1 1. These four equations were recast in the following matrix form Ax= b* 


fu e n e k 


0 -m e Ue( Ug.U^ Q ) 


ng£ ngirijUj 


Ui/ne 


0 kTe + meUeUi-meU^ 


where b is: 
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v ea n e k(T e -T a ) - 3 v e in e k(T e -Tj) - (£; + 2 ^T e + ^ rn e (U e 4-UgQ)) w - 

- mjUiw 
w/rie 

ne g U e 0 B r - n e m e Vea(U e -U a ) - n e meV e i(U e -Ui) - meU e w 



This system was then solved by recasting the matrix equation into a tridiagonal form and back 
substitutiing. This type of marching solution required initial conditions for each of the differential 
equation. The initial electron number density, electron temperature, electric potential and, atom and ion 
velocity were given. The values used for the initial conditions and the assumed values are shown in 
Table I. Due to numerical instabilities, the solution was restricted to positive temperature gradients in the 
initial portion of the discharge chamber even though negative gradients may exist. The use of a more 
sophisticated higher order solution technique might relax this constraint. 


Table I: Initial Conditions and Assumed Values For the Simplified Numerical Model 


Initial Conditions 

Assumed Values 

n e = lxlO 16 nr 3 

Ti= 1000K 

m = 4.5 mg/s 

T e = 5eV 

T fl = 1000K 

A = 0.0054 m2 

Ui = 500 m/s 

U a = 500 m/s 

$ rmax ^ 1 30 GHUSS 

<? = 0 Volts 

J = 4.5 Amperes 

UeQmax - lxlO 16 m/s 


Results and Discussion 

Analytical solutions were generated for a 3.6 cm long computational domain using the experimental 
collision cross sections. The length of the discharge chamber was 2 cm. These dimensions correspond 
to the geometry of the SPT-100. The radial magnetic field was given the Gaussian axial distribution 
shown in Figures 1 and 3. The azimuthal electron velocity was also assumed to have this functional 
dependance. This was deemed appropriate because the radial magnetic field is the primary mechanism for 
imparting azimuthal electron motion. The distribution was an analytic estimate of the actual magnetic field 
distribution. The peak value of the azimuthal electron velocity was calculated based on measured values 
of electron number density and azimuthal electron current[9,10]. The resulting potential distribution with 
the cathode potential set to zero, normalized assumed magnetic field distribution, electron number density 
and electron temperature are shown in Figure 3. 
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Figure 3: Spatial Distribution of the Calculated Potential, Electron Number Density and Electron Temperature and the input 
Radial Magnetic Field Distribution. 
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The potential dropped slowly until the magnetic field strength was sufficiently intense to significantly 
imped axial electron motion. In the region from z/L = 0.2 to 0.8 the potential dropped rapidly. Closer to 
the cathode plane the potential asymptoted to its final value. The total potential drop was approximately 
250 Volts which is comparable with the 300 Volt discharge voltage considered nominal for SPT-100 
operation. The electrode phenomena which were not included in this model may have contributed to this 
discrepancy. The maximum calculated axial ion velocities were 12000 m/s, which corresponds to a 
specific impulse of 1220 seconds, slightly below the nominal measured specific impulse of 1500 
seconds. 

The electron/ion number density increased gradually from the initial value of 10 16 nr 3 throughout the 
discharge until z/L = 0.8, beyond which it was nearly constant. The initial value corresponded to a 
ionization fraction of 0.001 or 0.1% at the anode. This value was chosen arbitrarily and a sensitivity 
analysis on this initial condition was not performed during this preliminary study. The peak value was 
approximately 3xl0 17 nr 3 . In terms of the percent of ionization, the plasma was essentially fully 
ionized by the exit, however, the ionization fraction didn't reach 0.90 until z/L =0.9. This indicated that 
ionization is occurring throughout the entire region of the discharge. The creation of the majority of 
electron/ion pairs at axial locations ranging from an z/L of 0.2 to 0.8 correspond to a location in the 
potential distribution where the majority of the potential drop occurred. This suggests that Joule heating 
of the electrons in the regions of high electric field strength, caused by significant radial magnetic field 
strengths, gave rise to effective ionization of neutral xenon in this region of the discharge which were 
subsequently accelerated away. 

The electron temperature was relatively constant in the initial portion of the discharge chamber where it 
was forced to have a positive gradient. Peak values of approximately 51 eV were obtained between the 
exit of the discharge chamber and the cathode plane. The calculated ionization rates had a measurable 
effect on this maximum value. When the rates generated using Drawin's empirical ionization cross 
sections were used a maximum electron temperature of 60 eV was calculated. Additionally, the 
ionization fraction at the exit dropped to 0.95 and the potential drop increased by 40 Volts to 290 Volts. 
Because of the excellent agreement between the various measurements of ionization cross sections these 
values were thought to be accurate. Consequently, the empirical values of ionization cross sections 
presented by Drawin, which had initially been adapted for convenience, were not used. 




Figure 4: Calculated Electron and Ion Current as a Function of Axial Position. Electron and Ion Currents at the Cathode 
Plane Showing Neutralization. 

Due to the reduced axial electron plasma conductivity the electrons were not the only significant current 
earners. As can be seen in Figure 4, the ions also carried a significant portion of the current. Near the 
anode the electrons essentially earned all the current. The amount of ion current increases and equals the 
electron current at an z/L of 0.7. In the remainder of the discharge the ions carried the majority of the 
current. However, the maximum ion current was restricted by the assumption of only singly ionized 
ions. For the supplied mass flow rate of 4.5 mg/s this corresponds to a maximum ion current of 3.3 
Amperes. What this implied was that the cathode supplied 1.2 Amperes to the discharge and the 
remaining 3.3 Amperes needed to neutralize the exiting ion beam as depicted in Figure 4. There is, 
however, the possibility during actual operation of ion beam currents in excess of the discharge current 
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due to the presence of multiply charged ions. This has previously been observed[ll]. In this instance 
the exiting ion beam would be neutralized by the electron cuiTent from the cathode, equal to the discharge 
current, and additional electrons from the discharge. Before these two sources of electrons, (/.<?. from the 
discharge and from the cathode), had experienced sufficient thermalizing collisions a non-equilibrium 
electron energy distribution function would exist consisting of the sum of two different equilibrium 
electron energy distribution functions. This is offered as a potential explanation of the non-equilibrium 
distribution experimentally observed [12,13]. If the electron energy distribution function was not 
Maxwellian within the discharge chamber, the rates used in the simplified numerical model developed 
during this investigation were invalid and these results do not apply. 



z/L 


Figure 5: Comparison of normalized Potential, 
Electron Temperature, and Electron Number 
Density with Fife[3] and Bishaev and Kim[9] 


Calculated solutions were compared to the those of Fife[3] 
and the measurements taken by Bishaev and Kim[9]. These 
comparisons are shown in Figure 5. In each of these 
investigations the discharge chamber was approximately 5 
cm long. In the SPT-100 the discharge chamber was 
approximately 2 cm long, therefore direct comparisons with 
these previous investigations required some form of 
normalization. The length of the discharge chamber was 
chosen. This distance was used rather than the anode to 
cathode distance since cathode position was not considered 
as critical. Furthermore, because the total potential drop 
was 200 Volts in the case of Fife and Bishaev and Kim the 
values were normalized by the maximum potential value. 
Additionally, their data were presented in two dimensions. 
Comparisons were made with their annular centerline 
values. The normalized potential distributions were similar 
with the majority of the potential drop taking place in a 
region near the exit of the discharge chamber and very little 
potential drop in the first half of the discharge chamber. 
Slight differences in the distribution may be attributed to the 
normalization scheme or differences in the thrusters being 
considered. 

The calculated values of electron temperature were 
considerably higher than those calculated by Fife and 
measured by Bishaev and Kim. The overall magnitudes 
were considerably higher. Both these results and the 
calculated results of Fife show a peak in electron 
temperature at or near the exit of the discharge chamber (z/L 
= 0.56). The measurements of Bishaev and Kim showed 
an electron temperature distributions with a peak midway 
between the discharge chamber exit and the cathode plane. 
Higher calculated exit electron maybe the result, of not 
including an expansion downstream of the constant area 
discharge chamber in the one-dimensional model presented 
here. 

The electron number densities are also compared. The 
calculated values and measurements between Bishaev and 
Kim show similar characteristics: peak values of 
approximately 3 to 5xl0 17 nr 3 between the discharge 
chamber exit and the cathode. The calculated values of Fife 
showed a substantially higher peak value occurring 
upstream of the midpoint of the discharge chamber. 
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Concluding Remarks 

A simplified numerical model of the plasma discharge within the SPT-100 stationary plasma thruster was 
presented. Seven coupled algebraic and differential equations, which included ion continuity, global 
continuity, current continuity, ion momentum, electron momentum and electron energy were shown to 
form a minimum subset of equations which retained all the essential physics, providing azimuthal electron 
velocities were included as an input parameter. Based on a consideration of the SPT-100 geometry with a 
4.5 mg/s anode flow, 4.5 Amperes of discharge current, and a peak radial magnetic field strength of 130 
Gauss, a quasi-one dimensional description used provided a solution set. This solution indicated peak 
plasma densities of approximately 3xl0 17 nr 3 between the discharge chamber exit and the cathode. Peak 
electron temperatures of 51 eV were calculated. This result was sensitive to the collision cross sections 
for ionization of xenon by electron impact. Ionization fractions of 1 were calculated near the cathode 
plane. The calculated potential distribution showed the majority of drop in the region of peak radial 
magnetic field strength, with a magnitude of 250 Volts. Calculated plasma densities and potential 
distribution compare favorably with previous measurements. Calculated electron temperatures are 
substantially higher than previously reported. 
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